%% This code runs Stage III of the method.
% assuming a linear cost trend
close all; clear;
clc;

% define saved files from running this code.
estFile_save = 'est2_output3_linearC.mat';
fig1_save = 'fig4c_supply_linearC.png';

%% load data and estimates from Stage I & II
load 'ATE_data_age10.mat'
load 'est1_output_2015.mat' % output year: 1950-2015

Tend = 2015;
ttmax = length(E);
H = (1e+3)*ATE_catch(year<=Tend); % harvest (kg)
Price = R16price(year<=Tend); % price ($/kg)

tt = year(1:end-1) - year(1);

%% Estimate coefficients in the dynamic Open Access difference eq
fprintf('DeltaE = gamma*(PH) + (gamma*(c0+c1*tt))(-E)\n');
% c(t) = c0+c1*t, linear cost trend

% generate variables for regression y = X*b
y = E(2:end)-E(1:end-1); % DeltaE
z = Price(1:end-1).*H(1:end-1); % Price*Catch
x0 = -E(1:end-1);
x1 = -tt.*E(1:end-1);

X = [z x0 x1];
fprintf('y = DeltaE, x1 = PH, x2 = (-E), x1 = t*(-E)\n');
result = fitlm(X,y,'Intercept',false);
disp(result)
[n, p] = size(X);
AIC = n*log(result.SSE/n)+2*p;

b3 = result.Coefficients.Estimate; % point estimates
se3 = result.Coefficients.SE; % std. err
gamma = b3(1);

%% Monte Carlo
M1 = 1000; % number of random draws
M2 = 500; % number of MC iterations

rng(54321) % for replication

% random sampling of gamma from Gamma dist.
pd1 = makedist('Gamma','a',(b3(1)/se3(1))^2, 'b',se3(1)^2/b3(1));
beta1_mc = random(pd1,M1,M2);
% random sampling of gamma*c0 from Gamma dist.
pd2 = makedist('Gamma','a',(b3(2)/se3(2))^2, 'b',se3(2)^2/b3(2));
beta2_mc = random(pd2,M1,M2);
% random sampling of gamma*c1 from Gamma dist.
pd3 = makedist('Gamma','a',(b3(3)/se3(3))^2, 'b',se3(3)^2/b3(3));
beta3_mc = random(pd3,M1,M2);
% sample mean of cost parameter
c0_sample = mean(beta2_mc./beta1_mc);
c1_sample = mean(beta3_mc./beta1_mc);

c0 = mean(c0_sample); se_c0 = std(c0_sample);
c1 = mean(c1_sample); se_c1 = std(c1_sample);

%% save cost and gamma estimates
save(estFile_save,'gamma','c0','se_c0','c1','se_c1')

%% Simulate backward-bending supply with linear c(t)
Rscale = R1_scale;

fig1 = figure('Position',[0 0 900 500]);
hold on
for t = 1:68 % simulate for each year 1950 to 2017
    cost = c0+c1*(t-1); % linear cost trend (mean estimate)
    [H_sim,P_sim,Hmax,backP] = sim_supply(Rscale,q0,Theta,cost);
    
    plot((1e-3)*H_sim,P_sim,'Color',0.8*[1 1 1],'LineWidth',1);
end
% (i) year = 1950
cost = c0;
[H_sim,P_sim,Hmax1,backP1] = sim_supply(Rscale,q0,Theta,cost);
plot((1e-3)*H_sim,P_sim,'Color',0.4*[1 1 1],'LineWidth',4);

% (ii) year = 2000
cost = c0 + c1*(2000-1950);
[H_sim,P_sim,Hmax2,backP2] = sim_supply(Rscale,q0,Theta,cost);
plot((1e-3)*H_sim,P_sim,'Color',0.2*[1 1 1],'LineWidth',3,'LineStyle',':');

% (iii) year = 2017
cost = c0 + c1*(2017-1950);
[H_sim,P_sim,Hmax3,backP3] = sim_supply(Rscale,q0,Theta,cost);
plot((1e-3)*H_sim,P_sim,'k-.','LineWidth',2.5);

xlabel 'Catch Quantity (metric tons)'
ylabel 'Price ($/kg)'
ylim([0 100])
xlim([0 Inf])
title '(c) simulated supply with linear cost trend'

% gtext(sprintf('1950: \n$%.1f/kg',backP1),'FontSize',12)
% gtext(sprintf('2000: \n$%.1f/kg',backP2),'FontSize',12)
% gtext(sprintf('2017: \n$%.1f/kg',backP3),'FontSize',12)
hold off

saveas(fig1,fig1_save)
